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ABSTRACT 

We present an optimized variant of the halo model, designed to produce accurate matter power 
spectra well into the non-linear regime for a wide range of cosmological models. To do this, 
we introduce physically motivated free parameters into the halo-model formalism and fit these 
to data from high-resolution /V-body simulations. For a variety of ACDM and vvCDM models 
the halo-model power is accurate to ~ 5 per cent for k < 10/zMpc 1 and z < 2. An advantage 
of our new halo model is that it can be adapted to account for the effects of baryonic feedback 
on the power spectrum. We demonstrate this by fitting the halo model to power spectra from 
the OWLS hydrodynamical simulation suite via parameters that govern halo internal structure. 
We are able to fit all feedback models investigated at the 5 per cent level using only two free 
parameters, and we place limits on the range of these halo parameters for feedback models in¬ 
vestigated by the OWLS simulations. Accurate predictions to high-k are vital for weak lensing 
surveys, and these halo parameters could be considered nuisance parameters to marginalize 
over in future analyses to mitigate uncertainty regarding the details of feedback. Finally, we in¬ 
vestigate how lensing observables predicted by our model compare to those from simulations 
and from HALOFIT for a range of k-cuts and feedback models and quantify the angular scales 
at which these effects become important. Code to calculate power spectra from the model 
presented in this paper can be found at https : //github. com/alexander-mead/hmcode. 

Key words: gravitational lensing: weak - cosmology: theory - dark energy - large-scale 
structure of Universe. 


1 INTRODUCTION 

In the standard theory of cosmological structure formation, all 
large-scale structure in the Universe forms via the gravitational col¬ 
lapse of small amplitude initial seed fluctuations. This process re¬ 
sults in a non-linear network of haloes, filaments and voids that is 
comprised of dark matter and baryons. One of the goals of modern 
cosmology is to probe these density fluctuations in the late Uni¬ 
verse and to use them to constrain models of the cosmos. In the 
early Universe, or at very large scales today, the density fluctua¬ 
tions are small in magnitude and can be analysed using linear per¬ 
turbation theory, which can be calculated precisely as a function 
of cosmological parameters - including both baryonic and dark- 
matter components (e.g., Seljak & Zaldarriaga 1996; Lewis et al. 
2000; Bias et al. 2011). As structures evolve in the later Universe 
they grow and become non-linear. Various perturbative schemes 
have been developed in cosmology to analyse these fluctuations 
(see reviews by Bernardeau et al. 2002; McQuinn & White 2015), 
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which give insight into the onset of non-linear structure formation. 
However, the most successful cosmological probes to date focus 
on the regime of linear perturbations, for example baryon acous¬ 
tic oscillations (BAO; e.g., Padmanabhan et al. 2012) or the cosmic 
microwave background (e.g., Planck Collaboration XIII 2015). In 
future surveys, the Universe will be mapped in finer detail and in 
principle it will be possible to extract a great deal of information 
from non-linear perturbations. 

Unfortunately, perturbative schemes fail as larger non- 
linearities develop, due to the inability of perturbation theory to 
model matter shell-crossing (see McQuinn & White 2015 for a ID 
discussion). Bound structures in the Universe today consist of mat¬ 
ter that has undergone many crossings and can represent large 
departures from the mean density. Currently, these extreme non- 
linearities can only be accurately modelled by running large cos¬ 
mological simulations, which commonly assume collisionless mat¬ 
ter, so-called dark-matter-only simulations. However, even accurate 
pure dark-matter simulations are computationally expensive and 
prohibit the wide space of possible cosmological parameters to be 
explored quickly. Furthermore, it can also be difficult to understand 
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which physical processes are at work in yielding a given simulation 
output. Thus, an analytic model for the evolution of structure can 
be invaluable, both in terms of speed and of insight. In this pa¬ 
per, we use the halo model (Peacock & Smith 2000; Seljak 2000; 
Cooray & Sheth 2002). This has become established as an impor¬ 
tant tool for explanation, but it is known not to provide the accuracy 
required for the interpretation of current data. 

A key measure of scale-dependent inhomogeneity, which can 
be calculated via perturbation theory or from the halo model, but 
also measured in simulations, is the power spectrum of the density 
field. Based on analytical insights, calibrated with A-body sim¬ 
ulations, various approximate formulae for the non-linear power 
spectrum have been generated. The most widely used of these to 
date has been the HALOFIT method of Smith et al. (2003), revised 
by Takahashi et al. (2012), which uses ideas from the halo model. 
This fitting formula has been expanded by various authors; to in¬ 
clude massive neutrinos (Bird et al. 2012) and /(/?) modified grav¬ 
ity models (Zhao 2014). 

A different approach that is used to make predictions for 
the non-linear spectrum is that of the ‘emulator’ code based on 
the ‘Coyote Universe" suite of simulations (Heitmann et al. 2009; 
Heitmann et al. 2010; Lawrence et al. 2010; Heitmann et al. 2014): 
the COSMIC EMU. Sets of high resolution simulations are run at key 
points in cosmological parameter space (so-called ‘nodes’) so as to 
cover the space evenly. The emulator then interpolates between the 
measured power spectra as a function of cosmology, yielding pre¬ 
dictions for any set of parameters within the space. Heitmann et al. 
(2014) show that their emulator produces the power spectrum to 
an accuracy of 1 per cent for k < IMpc -1 and 5 per cent for 
k < 10Mpc~* and it covers a small, but interesting, range of cos¬ 
mological parameter space for flat universes and dark energies with 
constant equation-of-state. However, COSMIC EMU makes no pre¬ 
dictions for k > lOfrMpc -1 or z > 4 and it would be useful to ex¬ 
pand the range of cosmological parameters that it currently encom¬ 
passes. 

Even if one were in possession of a perfect model of non¬ 
linear gravitational clustering it is difficult to compare this theory 
directly with the total matter field in the Universe. Instead one 
typically views the matter field via its gravitational light-bending 
effect, specifically in how bundles of light from a distant galaxy 
are sheared as they pass density perturbations between the galaxy 
and a telescope. This gravitational shearing induces a correla¬ 
tion in apparent galaxy shapes on the sky, as light from galax¬ 
ies that are close on the sky is distorted coherently. This ‘weak’ 
gravitational lensing (see the review by Kilbinger 2015) has been 
used to place constraints on cosmological models, for example 
COSMOS (Schrabback et al. 2010) or CFHTLenS (Heymans et al. 
2013; Kilbinger et al. 2013) and there are many forthcoming sur¬ 
veys designed to measure this effect in finer detail. Weak lensing 
measures a projected version of the total matter distribution in the 
Universe, in which the same shear correlations at a given angle can 
be caused by a smaller scale density fluctuation close to the ob¬ 
server, or a larger-scale fluctuation further away. This mixing of 
scales means that theoretical predictions for 2D weak-lensing ob¬ 
servables require predictions for the clustering of the full 3D mat¬ 
ter distribution over a wide range of scales and redshifts - although 
see the 3D lensing of Heavens (2003) and Kitching et al. (2014) for 
ways to avoid this. Current lensing analyses that work at the level of 
shear correlations either employ HALOFIT (e.g., Schrabback et al. 
2010; Heymans et al. 2013; Kilbinger et al. 2013), emulator-based 
strategies (Liu et al. 2015; Petri et al. 2015) or the technique of sim¬ 
ulation rescaling (Angulo & White 2010; Angulo & Hilbert 2015) 


to provide predictions for the matter spectrum for the required 
scales; these predictions are then converted to predictions for shear 
observables. 

In this paper, we take a completely different approach to mod¬ 
elling the full non-linear spectrum; we present an optimised variant 
of the halo model that is able to predict the matter power spectrum 
accurately to wave-numbers of interest for current and future lens¬ 
ing surveys (k ~ lO/iMpc -1 ). Our first goal is to provide accurate 
halo-model fits to dark-matter only simulations across a range of 
cosmological parameters; our approach is to identify parameters in 
the halo model that can be made to vary in a physically-motivated 
way and then to fit these to high-resolution simulated power spec¬ 
tra from the emulator presented in Heitmann et al. (2014). This 
approach is distinct from that of HALOFIT (Smith et al. 2003; 
Takahashi et al. 2012), which is an empirical fitting formula mo¬ 
tivated by the principles of the halo model, but one which does 
not use the halo model directly. Although we focus on weak lens¬ 
ing observables we stress that our optimised version of the halo 
model is general, and useful for any cosmological analysis that 
currently uses HALOFIT. Our approach is also distinct from recent 
work aimed at improving the halo model by Mohammed & Seljak 
(2014) and Seljak & Vlah (2015). These authors do not use the full 
apparatus of the halo model in order to provide accurate matches 
to the power spectrum, instead opting to employ a combination of 
perturbation theory and a series expansion, aimed primarily at an 
accurate modelling of the quasi-linear regime. 

An additional source of uncertainty is the impact of bary- 
onic feedback on the total matter distribution in the Universe. 
Most treatments of non-linear evolution ignore any interactions 
other than gravitational, except in the initial conditions. Work 
with perturbation theory at large scales (e.g., Shoji & Komatsu 
2009; Somogyi & Smith 2010) has shown that including the dis¬ 
tinct physics of dark matter and baryons offers small improve¬ 
ments (~ 0.5 per cent in power) compared to the approximation 
of treating ‘matter’ as being a single component. This has also 
been tested in simulations by Angulo et al. (2013), where it was 
shown that including distinct transfer functions for dark matter and 
baryons leads only to small differences in the eventual measured 
non-linear matter power spectrum - below one per cent at late 
times around the BAO scale. In contrast to these small effects at 
large scales, early semi-analytical treatments (White & Vale 2004; 
Zhan & Knox 2004) and recent work using hydrodynamical sim¬ 
ulations together with prescriptions for feedback (e.g., Jing et al. 
2006; Rudd et al. 2008; Schaye et al. 2010; Martizzi et al. 2014) 
have shown that the redistribution of matter caused by non-linear 
processes such as gas cooling, active-galaxy feedback and super¬ 
nova explosions can have a large impact on the total mass distribu¬ 
tion, but details regarding the magnitude of feedback are uncertain. 
We show that we are able to accurately capture the effects of a vari¬ 
ety of feedback recipes using our optimised halo model by varying 
only two parameters that govern halo internal structure. We con¬ 
strain these parameters and suggest limits for these that may form 
a prior. In forthcoming weak-lensing analyses these could either by 
constrained, to learn about feedback, or marginalized over to pro¬ 
vide unbiased cosmological constraints. 

This paper is structured as follows: In section 2 we first discuss 
our conventions for measures of inhomogeneity, then go on to dis¬ 
cuss the halo model, including the specifics of the particular model 
that we employ. Next, in section 3 we compare the original halo 
model to accurate A-body simulations to display its shortcomings 
and then discuss our modifications to the model. In section 4 we 
present our main result; the power spectra of the modified model 
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which best fits the emulator of Heitmann et al. (2014). We address 
baryonic feedback in section 5 and demonstrate that our approach 
can be easily adapted to account for feedback in a variety of mod¬ 
els. Finally in section 6 we show how our predictions for the matter 
spectrum translate into lensing observables and show differences 
between our approach, that of HALOFIT, and simulations. Our work 
is then summarized in section 7. Appendix A is dedicated to show¬ 
ing the response of the halo-model power spectrum to variations in 
cosmological parameters while Appendix B details the operation of 
our publicly available halo-model code. 


2 THE HALO MODEL 

2.1 Descriptions of inhomogeneity 

We will use the following notation: the matter density field, p, is 
given in terms of comoving position, x, and time, f, by 

P(x,f) = p(t)[l + S(x,t)] , (1) 

where p(f) is the mean matter density and 8(x,t) is the density 
fluctuation about the mean. We will be interested in the Fourier 
Space overdensity 4> which is defined via the Fourier convention 
of Peebles (1980) for a periodic volume V: 

* = £/ S(x)e~ ik ' x d 3 x ; (2) 


S(x) = £S k e' k *. 

k 


(3) 


The power spectrum of statistically isotropic density fluctuations 
depends only on k = |k| and is given by 

P{k) = (|4| 2 ) , (4) 


where the average is taken over modes with the same modulus but 
different orientations. We find it more convenient to use the dimen¬ 
sionless quantity A 2 : 

A 2 {k)m4nv(^j P(k), (5) 

which gives the fractional contribution to the variance per logarith¬ 
mic interval in k. If the field is filtered on some comoving scale R, 
the variance is 

g 2 (R)= [ A 2 (k)T 2 (kR)d\nk , (6) 

Jo 

where the window function is 


3 

T(x) = —j-fsinx — xcosx) , (7) 

x i 

in the case of smoothing with a spherical top-hat. 


2.2 Halo-model power spectra 

The halo model is an entirely analytic model for the non-linear 
matter distribution in the Universe that takes inspiration from the 
results of iV-body simulations. The basic idea here goes back to 
Neyman et al. (1953) but has been given a modern guise (Seljak 
2000; Peacock & Smith 2000), as reviewed in Cooray & Sheth 
(2002). The great power of the method is that it can encapsulate 
the clustering of galaxies by the choice of an appropriate halo- 
occupation number, specifying the number of galaxies as a function 


of halo mass. But the present application is simpler, since we are 
considering only the overall mass distribution. 

We now summarise the main features of the method, follow¬ 
ing most closely the presentation in Peacock & Smith (2000): The 
density field is described as a superposition of spherically symmet¬ 
ric haloes, with mass function and internal density structure that 
are accurately known as functions of cosmology from simulations. 
In the simplest case of randomly distributed spherical haloes the 
power spectrum has the form of shot noise, moderated by the den¬ 
sity profile of the haloes: 

A 2 m {k)=4n(^^j ^ J Q M 2 W 2 (k,M)F(M)dM . (8) 

Here the power spectrum is calculated as an integral over all halo 
masses, of mass M, where F(M) is the halo mass function (co¬ 
moving halo number density in d M) and W ( k,M ) is the normalized 
Fourier transform of the halo density profile: 

W(k,M) = — [ ) 4 jtr 2 p(r,M) dr , (9) 

M J o kr 

where r v is the halo virial radius. On large scales, haloes are not 
randomly distributed and displacements of haloes with respect to 
one another require us to consider a ‘two-halo’ term to the power. 
For the matter distribution, this is approximately the linear-theory 
power spectrum 

A 2 2H (k)=A 2 in (k). (10) 

An expression for the full halo-model power spectrum is then given 
by a simple sum of the terms 

A 2 (k) = A 2 ft + A. (11) 


2.3 Ingredients of the halo model 


To implement this model, we need to know the halo mass function 
and density profiles in order to evaluate the one-halo integral in 
equation (8). 

The halo density profile is commonly described via the 
Navarro, Frenk & White (1997; NFW) profile: 


P(r) 


_Pn_ 

(r/r s )(l+r/r s ) 2 ’ 


( 12 ) 


where r s is a scale radius that roughly separates the core of the halo 
from the outer portion and Pn is a normalization; in order to have 
a finite mass this profile must be truncated at the virial radius r v , 
within which the mean overdensity of the halo is A v . More recent 
work (Navarro et al. 2004) has shown that halo profiles can be bet¬ 
ter fitted with Einasto profiles, which differ from NFW near the 
halo centre. However, the halo-model power calculation (equation 
8) depends on a self-convolution of the profiles and this smears out 
details of the halo centre. Thus we prefer to use the simpler NFW 
fit. 

Simulated haloes need to be identified in a particle distribution 
and this is usually determined via a user-set overdensity threshold. 
Typically, a value of A v = 200 is taken, which is loosely based on 
predictions from the spherical collapse model in an = 1 uni¬ 
verse, although some authors use a value of A v that varies with cos¬ 
mological parameters in accordance with spherical model predic¬ 
tions (e.g., Bryan & Norman 1998). Once the over density thresh¬ 
old has been set and the halo mass measured the virial radius is no 
longer an independent parameter, and in order to conserve mass 


/ 3 M \ 1/3 
\4;rA v p ) 


(13) 
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Note that this means that in general the r v given by equation (13) 
will be different from the halo radius one may independently mea¬ 
sure from halo particles in a simulation. The normalization, Pn, is 
set by the requirement that the spherical integral of equation (12) 
gives the halo mass. The only free parameter in a fit to simula¬ 
tions is then r s , or equivalently the halo concentration c = r v /r s . 
An implication of this is that the value of c measured for simu¬ 
lated haloes depends on the halo definition used - particularly the 
A v criterion, the algorithm used to identify haloes and the scheme 
used for breaking up spurious haloes or unbinding particles (e.g., 
Knebe et al. 2011). 

Since the genesis of the NFW profile, a large number of re¬ 
lations between the concentration and mass of haloes have been 
developed. The general trend is that haloes of higher mass are 
less concentrated than those of lower mass, attributed to the fact 
that larger haloes formed in the more recent past and that the cen¬ 
tral density of a halo retains a memory of the cosmological den¬ 
sity at its formation time. The original c(M) relation proposed by 
Navarro et al. (1997) was shown to produce an incorrect redshift 
evolution by Bullock et al. (2001), who provided an updated rela¬ 
tion based around the concept of halo formation time. Around the 
same time a similar model by Eke et al. (2001) was introduced, 
which was intended to predict the correct c(M) relation in the 
case of models with the same background cosmological param¬ 
eters but different linear spectra, for example warm dark matter 
models compared to a cold dark matter (CDM) model. Lately focus 
has shifted to produce extremely accurate concentration-mass rela¬ 
tions for the standard ACDM cosmological model (e.g., Neto et al. 
2007; Gao et al. 2008) but these relations do not allow for gen¬ 
eral variations in cosmology. More recently, Prada et al. (2012) and 
Klypin et al. (2014) have suggested c(M) relations that are ‘univer¬ 
sal’, in that they do not depend on cosmology other than via the 
function cr(M) (equation 6). These relations predict that models 
with identical linear spectra should share a c(M ) relation, at odds 
with results from the concentration emulator of Kwan et al. (2013), 
which produces a different relation for models with identical linear 
spectra but different growth histories (e.g., a ACDM model com¬ 
pared to a wCDM model at z = 0 with identical tTs). 

We choose to use the relations of Bullock et al. (2001) be¬ 
cause it was derived by fitting to a wide variety of cosmologies and 
also because their haloes were defined with a cosmology-dependent 
overdensity criterion, and therefore naturally adapt to the changes 
that we plan to make to the halo model in section 3. The c(M) for¬ 
mula relates the concentration of a halo, identified at redshift z, to 
a formation redshift, Zf, via 


for £2 m = 1, with a very weak dependence on cosmology (see 
Eke, Cole & Frenk 1996 for flat models with A and Lacey & Cole 
1993 for matter-dominated open models). In Bullock et al. (2001) 
the parameters A = 4 and / = 0.01 were found from fitting the c(M) 
relation to halo profiles over a range of masses and cosmologies. 

For very massive haloes, equation (15) can assign a formation 
redshift that is less than the redshift under consideration, suggesting 
that the halo formed in the future. In our calculations we remedy 
this by setting c = A if Zf < z, although it makes very little practical 
difference to our power spectrum calculations. 

It was shown in Dolag et al. (2004) and Bartelmann et al. 
(2005) that the c(M ) relations proposed in Navarro et al. (1997), 
Bullock et al. (2001) or Eke et al. (2001) failed to reproduce the ex¬ 
act variations in concentration seen in models with identical linear- 
theory power spectra but different models of dark-energy. Differ¬ 
ences in concentration arise because haloes form at different times 
in these models, despite having matched linear theory at z = 0, and 
the exact form of this hysteresis was not being captured by existing 
relations (although the general trend is captured by Bullock et al. 
2001). Dolag et al. (2004) proposed a simple correction scheme 
that augments the ACDM concentration for a model by the ratio 
of asymptotic (z —> °°) growth factors of the dark-energy cosmol¬ 
ogy to the standard ACDM one: 


c de = ca 


gDE^-»°°) 

gA(z->°°) 


(16) 


and we implement this correction in our incarnation of the halo 
model. The effect of dark energy on halo concentrations can be 
seen at the level of the power spectrum in McDonald et al. (2006), 
in fig. 10 of Heitmann et al. (2014) and also in our Fig. Al. It can 
be seen at the level of measured halo concentrations using the c(M) 
emulator or Kwan et al. (2013). Because halo concentration affects 
small-scale power (equation 8), a corollary of this is that the full 
non-linear spectrum will be different at small scales in different 
dark-energy models, even if they share an identical linear spectrum. 
Any scheme in which the calculation of the non-linear power de¬ 
pends solely on the linear power will thus fail to capture this detail. 

The mass function of haloes (the fraction of haloes in the mass 
range M to M + AM) has been measured from simulations (e.g., 
Sheth & Tormen 1999; Jenkins et al. 2001) and has been shown 
to have a near-universal form, almost independent of cosmology, 
when expressed in terms of the variable 



(17) 


c(M,z)=A^, (14) 

1 +z 

where the parameter A is deduced by fitting to simulated haloes. 
The formation redshift is calculated by finding the redshift at which 
a fraction (/, also derived from simulated haloes) of the eventual 
halo mass has collapsed into objects, using the Press & Schechter 
(1974) theory: 

^ct(/M,z) = 5 c , (15) 

g(z) 

where g(z) is the linear-theory growth function normalized such 
that g(z = 0) = 1, cr“ is the variance of the linear density field fil¬ 
tered on the scale of a sphere containing a mass M (equation 6; M 
is the mass enclosed in a sphere with radius R in the homogeneous 
universe), and 5 C is the linear-theory collapse threshold. The value 
of 5 C is calculated from the spherical collapse model: 5 C ~ 1.686 


The mass function can be expressed as a universal function in /(v), 
which is related to F(M) that appears in equation (8) via 

- 1 -F(M)AM=f{v)Av . (18) 

This universality was predicted in an approach pioneered by 
Press & Schechter (1974) whereby the mass function was calcu¬ 
lated explicitly by considering what fraction of the density field, 
when smoothed on a given mass scale, is above the critical thresh¬ 
old for collapse (5 C ) at any given time. The expression that they 
calculated for the mass function is the Gaussian 



(19) 


but this is not a good fit to the mass function as measured in simu¬ 
lations; therefore we use the improved formula of Sheth & Tormen 
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Table 1. COSMIC EMU ranges for the six cosmological parameters that are 
allowed to vary. Note that ft), = £ljh 2 . £1^ and £2 m are the cosmological 
densities of baryons and all matter respectively, is spectral index of pri¬ 
mordial perturbations, h is the dimensionless Hubble parameter, w is the 
constant dark energy equation of state and as is the standard deviation of 
perturbations measured in the linear field at z — 0 when smoothed with a 
top-hat filter of radius 8/z~ 1 Mpc. Note that each model is kept flat, so that a 
change in ft)m at fixed h entails a change in £2 V , (the cosmological density in 
dark energy with equation of state w). The currently favoured Planck cos¬ 
mology (Planck Collaboration XIII 2015) lies close to the ‘node O’ values 
and the COSMIC EMU parameter space encompasses at least 5c deviations 
about this. Power spectra can be produced by the emulator between z — 0 
and 4 and for scales from k ~ 0.002/z to 10/iMpc -1 . 


Parameter 

Fiducial 

Minimum 

Maximum 

Node 0 

cob 

0.0225 

0.0215 

0.0235 

0.0224 

fflm 

0.1375 

0.120 

0.155 

0.1296 

«s 

0.95 

0.85 

1.05 

0.97 

h 

0.7 

0.55 

0.85 

0.72 

w 

-1 

-1.30 

-0.70 

-1 

og 

0.755 

0.616 

0.9 

0.8 


(1999), which was an empirical fit to simulations: 


/( v ) =A 


1 + 


(av 2 )P 


-av 2 /2 


( 20 ) 


where the parameters of the model are a = 0.707 and p = 0.3. A is 
constrained by the property that the integral of /(v) over all v must 
equal one, therefore A ce 0.2162. In Sheth & Tormen (1999) <5 C was 
taken to be 1.686, independent of the cosmology. 

Universality in the halo mass function is an approximation 
and the Sheth & Tormen (1999) mass function is only accurate at 
the ~ 20 per cent level (e.g., Warren et al. 2006; Reed et al. 2007; 
Lukic et al. 2007; Tinker et al. 2008; Courtin et al. 2011). How¬ 
ever, we do not attempt to use an updated non-universal mass func¬ 
tion prescription because we value the large parameter-space cover¬ 
age of Sheth & Tormen (1999) over a more accurate mass function 
tuned to only a small region of cosmological parameter space. 


3 OPTIMIZING THE HALO MODEL 

3.1 COSMIC EMU 

In this paper, we create an accurate halo model by fitting the halo 
model to data from high-resolution cosmological simulations. The 
simulation data we use comes from the COSMIC EMU, originally 
described in Heitmann et al. (2010) and updated in Heitmann et al. 
(2014). In these works, the authors ran a suite of high-resolution 
cosmological (V-body simulations and measured the power spec¬ 
trum in each. The simulations encompass a range of cosmological 
parameters (cOb, tOm, n s , h, w, erg) given in Table 1 and were con¬ 
structed to fill the parameter space in an ‘orthogonal Latin hyper¬ 
cube’ design to ensure that interpolating between models was as 
accurate as possible while still allowing only a small number (37) 
of simulations to be run. COSMIC EMU is the code released by the 
collaboration and allows the power to be produced at any point in 
their parameter cube from z = 0 to 4. The full spectra produced by 
COSMIC EMU are combinations of second-order Eulerian pertur¬ 
bation theory calculations at large scales and measured non-linear 
power at small scales. Extensive simulation resolution testing was 


conducted (Heitmann et al. 2010; Heitmann et al. 2014), and an ac¬ 
curacy of their simulations of 1 per cent in A 2 {k) to k = l/zMpc -1 
and 5 per cent to 10/iMpc^ 1 is quoted. Note that this is different 
from the accuracy of the interpolation scheme, which in the h -free 
version discussed in Heitmann et al. (2014) can be seen to be ~ 3 
per cent accurate to k = 1/zMpc -1 (degraded from the original 1 
per cent at k = l/zMpc~' when h is set free; see their figs 8 and 9) 
and ~ 5 per cent accurate to lO/zMpc -1 . 

Should we trust the accuracy stated for COSMIC EMU? Note 
that in answering this question it is important to specify whether 
one is comparing to the accuracy of the node simulations or of 
the emulator interpolation method. The nodes of COSMIC EMU 
are compared to simulations in Takahashi et al. (2012) where it is 
shown that a subset of the simulations of COSMIC EMU (MOOl to 
M009 in Lawrence et al. 2010) seem to have systematically ~ 2 per 
cent less power at k = 1/zMpc -1 when compared to the simulations 
presented in Takahashi et al. (2012). However, the authors do not 
present their own resolution tests, and resolution issues may plau¬ 
sibly be the origin of this small discrepancy. In Smith et al. (2014) 
the accuracy of power spectrum predictions was investigated as a 
function of ‘non-physical’ simulation parameters (such as force and 
mass resolution). The authors report that changing the PMGRID pa¬ 
rameter in GADGET-2 (Springel 2005; the code used for the COS¬ 
MIC EMU simulations) can have the surprisingly large effect of ~ 3 
per cent power differences at k = l/iMpc~*. This parameter was 
not investigated in the resolution testing of Heitmann et al. (2010). 
More recently Skillman et al. (2014) and Schneider et al. (2015) 
have shown disagreements around the 3 to 5 per cent level when 
comparing their own simulations to COSMIC EMU for the Planck 
cosmology (not one of the emulator nodes) for k < lO/zMpc -1 , 
but their claims are based on their own simulations of a single 
cosmological model and the error is within that quoted for the li¬ 
tres Heitmann et al. (2014) extended emulator. Given this discus¬ 
sion we err on the side of trusting the stated simulation accuracy 
in the COSMIC EMU papers with the caveat that the PMGRID claim 
of Smith et al. (2014) warrants further investigation and that addi¬ 
tional comparisons to different suites of simulations would be ben¬ 
eficial. 

We also note the existence of another power spectrum pre¬ 
diction code, PKANN (Agarwal et al. 2012; Agarwal et al. 2014), 
which uses neural networks to carry out interpolation between sim¬ 
ulation nodes. We choose not to use this because the simulations it 
was trained on are not as high resolution as those of COSMIC EMU 
(limited to k < 1/zMpc -1 ). In addition, the PKANN simulations con¬ 
tain some hydrodynamics and we wanted to initially focus only on 
dark matter. However, PKANN does allow the user to vary neutrino 
mass, which may be useful in further work. 

Although having accurate power spectra is extremely useful, 
three obvious limitations of the COSMIC EMU method are: How to 
extend the predictions to k> lO/zMpc -1 , or z > 4? How to extend 
beyond the emulator cosmological parameter space? And how to 
account for baryonic feedback? The first question arises in weak- 
lensing studies where predictions for lensing observables techni¬ 
cally require integrals of A 2 over all k (see section 6) but the main 
challenge in extending to smaller scales is to include the complica¬ 
tions of baryon feedback. However, modelling these scales is nec¬ 
essary if weak lensing is to achieve its claimed future precision. 
Additionally, in Monte Carlo Markov chain (MCMC) analyses, the 
chains will inevitably wish to explore outside the parameter range 
of COSMIC EMU and it is unclear how to proceed in this case. In this 
paper, we fit a variant of the halo model to data from the ‘nodes’ of 
COSMIC EMU, which are the exact locations within the cosmolog- 
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Figure 1. A comparison of the original halo model, described in section 2.2, 
to node 0 of COSMIC EMU (£2 m = 0.25, £2 b = 0.043, n s = 0.97, w = —1, 
(7* = 0.8, h = 0.72) at z 0 (upper; solid red), and z = 1 (lower; solid 
blue). The spectrum from COSMIC EMU is shown as black crosses. One can 
see that the halo-model power spectrum is qualitatively correct in shape, 
but under predicts the true power at the tens of per cent level at scales k > 
0.2/iMpc at both redshifts. This is exactly the scale at which the one-halo 
term comes to dominate over the linear power (short-dashed lines in the top 
panel). The long-dashed lines show a preview of the final fits we are able to 
produce in this work, which agree with COSMIC EMU at the 5 per cent level 
across all scales. In the top panel, our model seems to track the simulation 
almost exactly and is hard to distinguish. 


ical parameter space where the simulations were run. This has two 
advantages: The accuracy of the emulator is likely to be highest at 
the nodes, because there is no interpolation taking place. Secondly, 
in using the nodes we are taking advantage of the Latin hypercube 
deign of COSMIC EMU. Our resulting halo model fits can be used 
to extend the simulations to higher k, or higher z, in a physical way 
because they are motivated by theoretical arguments. Additionally, 
we show in section 5 that the halo model can be adapted to account 
for the influence of baryons on the matter power spectrum, by fit¬ 
ting parameters relating to halo internal structure to data from hy- 
drodynamical simulations. We also suggest that a successful fitting 
recipe can be used directly to explore models outside the COSMIC 
EMU parameter space. 

In Fig. 1 , we show a comparison of the power spectrum pre¬ 
dicted by our original incarnation of the halo model (NFW haloes; 
Bullock et al. 2001 concentrations; Sheth & Tormen 1999 mass 
function; A v = 200; <5 C = 1.686) to the power spectrum of COS¬ 
MIC EMU node 0, which is vanilla ACDM near the centre of the 
parameter space, at z = 0 and 1. It is immediately obvious that the 
halo-model prediction is qualitatively reasonable in form, but de¬ 
viates in detail from the simulations showing an underestimate of 
power of ~ 30 per cent for k > 0.5/iMpc -1 . There are several pos¬ 
sible reasons for the relatively poor performance of the halo model; 
Halo-finding algorithms tend only to assign half of the particles 


in a simulation into haloes (lenkins et al. 2001; More et al. 2011) 
so the non-linear distribution of half of the mass in the simulation 
is treated by the halo model via an extrapolation of the formula 
for the mass function. There are also clearly unvirialized objects in 
the quasi-linear regime that are not taken into account in our halo- 
model formalism, which also neglects halo substructure and as- 
phericity as well as non-linear material that may lie outside the halo 
virial radius. In addition, a scatter in any halo property at fixed mass 
will change the predicted halo-model power spectrum. For exam¬ 
ple, Cooray & Hu (2001) investigate a halo model with a scatter in 
c(M), which typically boosts the power, while Giocoli et al. (2010) 
also include the power due to halo substructure via a substructure 
mass function. How the measured power spectrum is altered under 
various assumptions can also be seen in recent simulation work by 
van Daalen & Schaye (2015) or Pace et al. (2015). 

Other problems are visible at large scales, where the halo 
model power can be seen to over-predict the simulations for k < 
0.1/7Mpc _1 . At these scales, the power is mildly non-linear and 
the two-halo term is in error, as well as the two- to one-halo tran¬ 
sition. Attempts to accurately model quasi-linear scales using the 
halo model have been made by Valageas & Nishimichi (2011), 
Mohammed & Seljak (2014) and Seljak & Vlah (2015), who use 
perturbation theory results as a two-halo term, and by Smith et al. 
(2007) who includes non-linear halo bias in the two-halo term. Ac¬ 
curate modelling of mildly non-linear power is an active field of 
research due to the importance of these scales for BAO measure¬ 
ments. 

3.2 Fitting a general halo model 

Rather than attempting to improve the halo model by adding miss¬ 
ing ingredients (e.g., Smith et al. 2007; Giocoli etal. 2010), thus 
making it more complicated, in this paper we take a more prag¬ 
matic approach: it is possible that part of the inaccuracy of the 
power spectrum calculation stems partly from incorrect parame¬ 
ter choices. The model contains quantities such as A v , which are 
round numbers motivated by analytic arguments. We may there¬ 
fore hope that improved results may be obtained by fitting the 
halo model to simulated power spectra using these quantities as 
physically-motivated free parameters. Our proposed changes rep¬ 
resent a prescription for producing effective haloes whose power 
spectrum mimics the true one, even if these haloes differ from 
those measured directly in simulations. The hope is that we can 
trade off inaccuracies in e.g., halo concentration against issues that 
are neglected in the standard halo model (asphericity, substructures, 
scatter in halo profiles), such that the two-point predictions are im¬ 
proved. 

Nevertheless, we wish to retain the large amount of tested the¬ 
oretical input that goes into the halo model. For example: changes 
in cosmological parameters alter the linear power spectrum, which 
in turn affects the mass function through the variance and the halo 
density profiles through the concentration and size relations. In ad¬ 
dition, the linear growth rate will change, which also affects the 
concentration relations directly as well as the amplitude of the 
linear power spectrum. Since all of these ingredients have been 
tested against simulations, there are grounds for hoping that a small 
amount of parameter readjustment may allow the halo model to 
produce robust predictions for the non-linear power spectrum that 
are of useful accuracy for a wide range of cosmological parame¬ 
ters. The dotted lines in Fig. 1 give a preview of how fruitful this 
approach is in providing an accurate model of the non-linear matter 
power spectrum. 
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3.2.1 Adapting the two-halo term 

The two-halo term governs power on large scales and is given in its 
original form in equation (10). Linear theory slightly over-predicts 
the matter power spectrum around the quasi-linear scale and does 
a particularly poor job of modelling damping of the BAO peaks at 
z = 0 , which are damped by the quasi-linear effect of small-scale 
displacements. Modelling of the minutiae of the damping of the 
BAO peaks is beyond the scope of this work, but we treat the damp¬ 
ing around these scales based on a model for the damping predicted 
from perturbation theory by Crocce et al. (2006), where 


A t in (k)^ e -^A^ n (k), (21) 

and <Ty is the 1 D linear-theory displacement variance given by 





u: 


k 2 


d k . 


( 22 ) 


The derivation of this expression assumes that the scales of interest 
are large compared to <J W , such that the damping factor cannot be 
trusted when k( 7 V is large. We found that the best fit to numerical 
data at this scale required an expression equal to equation ( 21 ) to 
quadratic order, but without the extreme high-Ac truncation: 


^m(k) = [l -/tanh 2 (tov/v 7 /)] > (23) 

where / is a free parameter in our fit. In the k( 7 V 3> 1 limit, equa¬ 
tion (23) reduces to A 2 H = (1 — /)A 2 in . 


3.2.2 Adapting the one-lialo term 

We add freedom to the canonical form of the one-halo term in equa¬ 
tion ( 8 ) in a number of ways: The first concerns the behaviour 
of the one-halo term at large scales, where the Universe tends to 
homogeneity faster than predicted by Poisson shot noise. At large 
scales, the one-halo term in equation ( 8 ) decays as A 2 ock 2 , whereas 
the linear power decays approximately °c k 4 , so it is inevitable that 
the one-halo term becomes greater than linear theory on very large 
scales, which is unphysical. This effect arises because haloes are 
treated as randomly placed in the standard halo-model formalism, 
when in fact they are clustered and distributed more smoothly than 
uniform random on very large scales. It has been suggested that a 
large-scale cut-off in the one-halo term can be physically explained 
as ‘halo exclusion' (Smith et al. 2007) an effect that arises because, 
by definition, haloes cannot exist within each other. This is not cap¬ 
tured by the standard halo-model power calculation because that 
calculation assumes that haloes are randomly placed, so that the 
probability of haloes overlapping is non-zero. Accounting for halo 
exclusion damps the halo-model power on large scales. Regardless 
of the exact details of exclusion, we modify the one-halo term so 
that it decays more rapidly than linear theory at large scales: 


A,i, - P 




■]Af 


H i 


(24) 


where Ay H is the same as in equation ( 8 ) and k * is a free parameter. 

Within the one-halo term, parameters that we allow to vary 
are the virialized overdensity of a halo, A v , defined in equa¬ 
tion (13), and the linear collapse threshold, 8 C , defined in equa¬ 
tion (17). Both of these parameters derive from the spherical model 
(e.g., p. 488 of Peacock 1999) and rely on a somewhat arbitrary 
definition of the exact time of halo collapse. The variation of 
A v can be predicted theoretically from the spherical model, and 
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Bryan & Norman (1998) provide a fitting formula 1 for a ACDM 
cosmology 


A v = 


18?r 2 + 82[Q m (z) — 1] — 39[Q m (z) — l ] 2 
£2m(z) 


(25) 


This suggests that A v increases as the universe deviates from £2 m = 

1 . 

In standard theory, 8 C ~ 1 .686 but we allow this number to be a 
free parameter in our fit to power spectrum data. Note that changing 
S c changes the relationship between v and the halo mass (equation 
17). This means that the ‘effective’ mass function we invoke to 
improve A 2 (k) predictions will not necessarily accurately represent 
the mass function that might be measured in simulations. 

Fitted halo relations, such as the mass function and mass- 
concentration relation, depend upon how haloes are defined when 
identified in simulations. Therefore, the variations of A v in our fit¬ 
ted halo model may not follow the simple theoretical variation in 
equation (25) exactly, but we assume that the trend of increased A v 
as the universe deviates from f 2 m = 1 will serve as a useful initial 
guide when we explore parameter space. In addition, for flat models 
with a single component of dark energy it is expected that A v would 
be a function of £ 2 m (z) only and this will be a useful principle in 
parameterizing fitting formulae. Increasing A v has the effect of in¬ 
creasing the internal density of haloes and thus decreases the virial 
radius of a halo of a fixed mass, thus increasing small-scale power. 
Increasing 8 C means the linear density field has to reach higher val¬ 
ues before collapse can occur (in the Press & Schechter 1974 ap¬ 
proach), the result of which is that the density field is dissected into 
more haloes of lower mass, which will reduce the amplitude of the 
shot-noise component of the one-halo term and thus reduce power. 

One further free parameter is 7], which we use to alter the halo 
window function via 


W{k,M) ->W(v T) fe,M) , (26) 

changing the halo profile in a mass dependent way but leaving 
V = 1 haloes unaltered and the individual halo masses unchanged. 
For rj > 0 higher mass (v > 1) haloes are puffed out while lower 
mass haloes are contracted, both at constant virial radius: 17 > 0 
decreases the power whereas r) < 0 increases it. This extra ingredi¬ 
ent was introduced to control the curvature of the power spectrum 
beyond k ~ 1 /tMpc^ 1 , where the filtering effect from the typical 
haloes has a major effect on the shape of the one-halo term. As we 
move to higher k values, the properties of lower mass haloes be¬ 
come increasingly important. It is difficult for the one-halo term to 
track to the smallest scales, and correcting this requires an empir¬ 
ical perturbation of the halo profiles. Additionally, we allow our¬ 
selves to vary the amplitude of the concentration-mass relation: A 
in equation (14). 


3.2.3 Full power 

A well known defect in the halo model is in the transition between 
the one- and two-halo terms, the so-called quasi-linear regime. In 
the standard halo model, the transition is modelled by a simple 
sum of the one- and two-halo terms (equation 11 ), but this is ob¬ 
viously deficient. At z = 0, this transition scale is approximately 
k = O.lh Mpc~* corresponding to physical scales of the order of 


1 Equation (25) differs slightly from that in Bryan & Norman (1998) be¬ 
cause we work with respect to the matter density, rather than critical den¬ 
sity. 
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Figure 2. Halo-model fits to all 37 nodes from COSMIC EMU at z = 0 (left) 0.5, 1 and 2 (right) are shown in the top row, while the bottom row shows the 
predictions from the Takahashi et al. (2012) revision of HALOFIT. The average fit is shown as the thick black line, whereas each individual thin coloured line 
shows a specific node. The range of the coloured lines in each panel give an idea of how the accuracy of the model varies with cosmological parameters. The 
two models are of comparable accuracy for these cosmological models shown, although the halo-model approach performs slightly better. The halo-model 
fitting formula performs worst at z = 2 where the high-k power becomes systematically inaccurate and at z = 0 where there is a spread in power at high k 
values and a systematic under prediction in power around k = 0.3/zMpc . However, it is mainly accurate at the 5 per cent level for all models and all scales 
shown. COSMIC EMU itself is claimed to be accurate at the 1 per cent level at k — 1/zMpc -1 and 5 per cent at lOfiMpc -1 , so our fit is mainly within this error 
around k ~ 10/iMpc -1 . 


Table 2. Halo-model parameter descriptions and values before and after fitting 


Parameter 

Description 

Original value 

Fitted value 

Equation in text 

A v 

Virialized halo overdensity 

200 

418 x £2~ 0352 (z) 

13 

4 

Linear collapse threshold 

1.686 

1.59 + 0.0314 In <r 8 (z) 

17 

0 

Halo bloating parameter 

0 

0.603 — 0.3(78 (z) 

26 

/ 

Linear spectrum transition damping factor 

0 

0.188 x o% 29 (z) 

23 

k t 

One-halo damping wavenumber 

0 

0.584 x cr v -1 (z) 

24 

A 

Minimum halo concentration 

4 

3.13 

14 

a 

Quasi-linear one- to two-halo term softening 

1 

2.93x1.77*" 

27 


tens of Mpc. On these scales, contributions to the density field will 
include, but are not limited to, large structure at the turn-around ra¬ 
dius, sheets, filaments and voids. It would be rather surprising if the 
complexity of non-linear evolution on these scales could be accu¬ 
rately modelled by a simple sum of crude one- and two-halo terms. 
In testing, we noted that the halo model performed most poorly 
around these transition scales and we address this problem by mod¬ 


elling the transition via 

A 2 (*) = [(A' 2 2 H )“ + (A' 1 2 H )“] 1 /“, (27) 

where a is the final parameter that we adjust to match simula¬ 
tions. Values of a < 1 soften the transition between the two terms 
whereas a > 1 sharpen it. The power at these scales is quite smooth, 
so fitting the transition via a is sufficient. 
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Table 3. Cosmological parameters inferred from various data analyses. In 
all cases, we quote the best fit with w = — 1 and flatness enforced. 


Cosmology 

fib 

iQ m 

n s 

h 

08 

WMAPl 

0.0457 

0.275 

0.969 

0.701 

0.810 

WMAP 9 

0.0473 

0.291 

0.969 

0.690 

0.826 

CFHTLenS 

0.0437 

0.255 

0.967 

0.717 

0.794 

Planck EE 

0.0487 

0.286 

0.973 

0.702 

0.796 

Planck All 

0.0492 

0.314 

0.965 

0.673 

0.831 


4 RESULTS 


We fit the parameters introduced in the previous section to data 
from all 37 nodes of COSMIC EMU at redshifts z = 0, 0.5, 1, 1.5 
and 2 with equal weight given to each redshift and node and k 
weighted equally in logarithmic space from 0.01 to 10/; Mpc~ 1 . We 
use a least squares method to characterize goodness-of-fit and use 
an MCMC-like approach to fit all parameters simultaneously. Our 
best-fitting parameters are given in Table 2 where there are a total 
of 12 parameters that are fitted to simulations, which can be com¬ 
pared with 34 for the Takahashi et al. (2012) version of HALOFIT. 
The cosmological dependences of each of our parameters was in¬ 
ferred by some experimentation. In Table 2, we see that a depends 
on n e ff, which is the effective spectral index of the linear power 
spectrum at the collapse scale, defined in Smith et al. (2003): 


3 + n eff = 


dlna 2 (R) 

dlnR 


< 7=1 


(28) 


However, our n e ff is slightly different from that in Smith et al. 
(2003) because we define u(R) using a top-hat filter, rather than 
a Gaussian. 

The accuracy of this model is demonstrated in the upper row 
of Fig. 2, which shows a ratio of the halo model to COSMIC EMU 
at z = 0, 0.5, 1 and 2. One can see that our fitted halo-model 
predictions are mainly accurate to within 5 per cent across all 
redshifts for the range of scales shown. We call this calibrated 
halo model HMCODE and refer to it thus throughout the remain¬ 
der of this work. We also tested our model at z = 3, a redshift 
to which it was not calibrated, and found that errors rarely ex¬ 
ceed 10 per cent. Takahashi et al. (2012) use the framework of the 
original HALOFIT of Smith et al. (2003), but obtain improved ac¬ 
curacy by fitting to modem simulation data with superior resolu¬ 
tion, extending to k = 30/;Mpc~’. The authors also focus their 
attention on models close to the current ACDM paradigm, rather 
than more general models (such as those with power-law spec¬ 
tra, or curved models). Takahashi et al. (2012) used simulations 
of 16 different cosmological models around the best fits from the 
Wilkinson Microwave Anisotropy Probe ( WMAP ) satellite (WMAP 
7 - Komatsu et al. 2011; WMAP 9 - Hinshaw et al. 2013) and in¬ 
clude models with tv^-1. One can see how well Takahashi et al. 
(2012) compare to COSMIC EMU in the lower row of Fig. 2 where 
HALOFIT can be seen to be comparable to our halo model but there 
is more high-/: spread at z = 0 and a systematic over-prediction 
of the power around k = l/;Mpc~* that worsens with increasing 
redshift. The stated accuracy of this version of HALOFIT is 5 per 
cent for k < 1/zMpc -1 and 10 per cent up to 10/zMpc , which 
is consistent with what is seen here. A similar plot for the orig¬ 
inal Smith et al. (2003) version of HALOFIT shows large under¬ 
predictions for k > 0.5/;Mpc~*. From this point onwards we only 
compare to the revised Takahashi et al. (2012) version of HALOFIT. 

In Fig. 3 we show how our model fares for cosmological pa- 



0.01 0.1 1 10 
k/(h Mpc' 1 ) 


Figure 3. A comparison of the power spectrum at z = 0.5 of HMCODE and 
HALOFIT to that of COSMIC EMU for several commonly used cosmological 
models (see Table 3) that derive from recent data sets. The error for each 
model is very similar because the cosmological models are all relatively 
similar. The HMCODE error rarely exceeds 2 per cent, with the exception 
being around the BAO peak, which stems from our not modelling the non¬ 
linear damping of the BAO. The HALOFIT error rises to around 4 per cent 
for k > 1/zMpc for all models. 


rameters derived from recent data sets (see Table 3). Once again 
we compare to COSMIC EMU and show results for both our cali¬ 
brated halo model and for the Takahashi et al. (2012) HALOFIT at 
z = 0.5. One can see that the error from the halo-model approach 
rarely exceeds 2 per cent for k < lO/iMpc -1 for these cosmologies, 
with the worst error being an over prediction of the amplitude of 
the BAO peaks around k = 0.2/; Mpc . This arises because we did 
not attempt to model the exact non-linear damping of this feature in 
the power spectrum, and so our prediction here is very close to the 
undamped linear prediction. That our errors are better here than for 
the more general models shown in Fig. 2 is because these models 
all lie close to the centre of the COSMIC EMU parameter space (see 
Table 1). The Takahashi et al. (2012) HALOFIT model works bet¬ 
ter at BAO scales, but over-predicts the power at k > 0.5/;Mpc~* 
systematically at around the 4 per cent level. 

The model presented here performs similar to, but slightly bet¬ 
ter than, the Takahashi et al. (2012) version of HALOFIT and has 
several advantages. Foremost, because we retain the apparatus of 
the halo model in our calculation, it means we can produce A 2 ( k ) 
to arbitrarily high k in a physically motivated way. Even though 
such extreme scales receive a small weight in lensing, they can 
be important if the modelling is badly wrong in this regime (e.g., 
Hamois-Deraps et al. 2015). A polynomial-based fitting formula 
such as HALOFIT risks generating pathological results when mov¬ 
ing beyond the regime constrained by simulations and it is not at 
all obvious how to extend COSMIC EMU. In Fig. 4, we show a com¬ 
parison of the power spectrum predicted out to k = lOO/zMpc -1 
with different models, simply to illustrate the range of behaviour at 
k > lO/iMpc -1 . Given that no simulations exist that could claim to 
accurately predict the matter power spectrum to k = lOO/zMpc -1 
at z = 0 we cannot make any quantitative statements about the ac¬ 
curacy of either model at these extreme scales, although both per¬ 
form comparably. The grey shaded region in Fig. 4 delimits these 
extreme scales and it is interesting to note that the maximum devi¬ 
ation between our model and HALOFIT is only ~ 10 per cent. 
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Figure 4. The predictions of HMCODE compared to the two commonly used 
HALOFIT schemes up to k = 100/zMpc -1 for a Planck cosmology at z = 0.5. 
The upper panel shows A 2 (k), while the lower panel shows the ratio of 
HMCODE to Takahashi et al. (2012); we cannot show a comparison with 
COSMIC EMU because it makes no predictions beyond k — 10/zMpc -1 . The 
range outside the bounds of COSMIC EMU is marked by the grey region. 
The HMCODE and Takahashi et al. (2012) models make predictions within 
5 per cent to k — 10/iMpc -1 but this discrepancy increases to 10 per cent 
for k < 100/iMpc 1 . Note that the power at these small scales is certainly 
strongly influenced by baryonic physics. The general level of agreement 
between these two models for a range of cosmologies can be inferred from 
Fig. 2. 

Recently an approach related to ours has been pursued by 
Mohammed & Seljak (2014) and Seljak & Vlah (2015); these au¬ 
thors use a combination of perturbation theory for a two-halo term 
and a polynomial series expansion for a one-halo term, constrained 
to contain only even powers of k by theoretical considerations. 
These authors obtain remarkable fits to the matter power spectrum 
and correlation function, but at the cost of fitting each term in the 
one-halo series expansion to simulations up to k ~ l/iMpc -1 . In 
Mohammed & Seljak (2014), it was shown that fits accurate at the 
2 per cent level were possible to the COSMIC EMU nodes up to 
k = 0.3/;Mpc~* but their fit degrades at smaller scales; it was tested 
out to 0.7/iMpc _1 and it is not obvious how to extend predictions to 
smaller scales. Our approach utilizes the full apparatus of the halo 
model and can therefore be extrapolated with a degree of robust¬ 
ness to the smaller scales that are required to make predictions for 
weak-lensing observables, which we show in section 6. We are also 
in a position to be able to model the effects of baryonic feedback 
in a physically motivated way, and we now turn the attention of the 
reader to this subject. 


5 BARYONIC PHYSICS 

Whilst 5 per cent accuracy across a range of cosmologies and scales 
is an important achievement of this work, it is clear that baryonic 
processes can have a much larger impact on the non-linear mass 
distribution than incorrect modelling of the dark-matter only spec¬ 
trum (e.g., White & Vale 2004; Zhan & Knox 2004; Zentner et al. 
2008; Casarini et al. 2011; Zentner et al. 2013; Semboloni et al. 
2013; Eifleretal. 2014; Hamois-Deraps et al. 2015). Baryonic 
physics is not normally accounted for in numerical simulations and 
baryons can undergo processes such as radiative cooling where they 


collect in sufficient density. The gas then contracts, which alters 
the dark matter distribution via gravitational interactions, and so 
the total matter distribution is altered because neither the baryons 
nor dark matter are where they would be if only gravitational inter¬ 
actions were considered (e.g., Jing et al. 2006; Duffy et al. 2010). 
Alternatively, supernova explosions or energy released by active 
galactic nuclei (AGN) can heat gas, which can then expand outside 
of the virial radius of haloes (Schaye et al. 2010; van Daalen et al. 
2011; Martizzi et al. 2014; van Daalen & Schaye 2015), and as a 
consequence the total matter distribution within a halo can be al¬ 
tered significantly. AGN feedback can reduce the baryon fraction 
in the centres of haloes by a factor of 2 in the most extreme models 
(Duffy et al. 2010). 

An advantage of the approach advocated in this work is that 
one might attempt to capture the influence of baryons on the mat¬ 
ter power spectrum simply by varying parameters that control the 
internal structure of haloes. This is possible because we retain 
the theoretical halo-model apparatus in HMCODE. Physically, one 
can regard baryonic processes as altering the internal structure of 
haloes, while not affecting their positions or masses to the same 
degree (e.g., van Daalen et al. 2014). It has been demonstrated that 
the effect of baryons can be captured by altering the halo internal 
structure relations, using information that is measured in baryonic 
simulations (Zentner et al. 2008; Duffy et al. 2010; Zentner et al. 
2013; Semboloni et al. 2013; Mohammed et al. 2014). The general 
trend is that gas cooling increases the central density of haloes 
whereas violent feedback, such as that from AGN, decreases the 
concentration. How this translates into the matter power spectrum 
in simulations is considered in van Daalen et al. (2011) where it 
was shown that per cent level changes in A 2 (k) can arise at k = 
0.3/iMpc -1 as a result of strong AGN feedback. Semboloni et al. 
(2011), Eifler et al. (2014) and Mohammed et al. (2014) all showed 
that failing to account for feedback would strongly bias cosmo¬ 
logical constraints from the weak lensing Dark Energy Survey if 
the most extreme feedback scenarios apply to our Universe and 
constraints from Euclid would be severely biased for any feasi¬ 
ble feedback scenario. Duffy et al. (2010), Semboloni et al. (201 1) 
and Mohammed et al. (2014) also showed that the main effects of 
baryonic feedback could be captured using a halo-model prescrip¬ 
tion, considering how feedback would alter the internal structure of 
haloes. 

We use power spectra from the Overwhelmingly Large Simu¬ 
lations (OWLS; Schaye et al. 2010; spectra from van Daalen et al. 
2011) of a dark-matter (DMONLY) model; a model that has pre¬ 
scriptions for gas cooling, heating, star formation and evolution, 
chemical enrichment and supernovae feedback (REF); a model 
that is similar to REF but with the addition of active galactic 
nuclei (AGN) feedback (called AGN); and a model similar to 
REF but which additionally has a top-heavy stellar initial mass 
function and extra supernova energy in wind velocity (DBLIM— 
called DBLIMFV 1618 in van Daalen et al. 2011). It was shown in 
van Daalen et al. (201 1) that the difference in power between the 
DMONLY and AGN models is particularly large. 

We fit the power spectra from the OWLS simulations using our 
calibrated halo-model approach with a halo profile that is altered 
to reflect baryon bloating and gas cooling. Again, our new fitted 
halo profiles may not match those of simulated haloes exactly but 
our aim is to match the power spectrum accurately. However, we 
would expect the trends observed in the profiles of simulated haloes 
to be respected by any modification to halo profiles in HMCODE. 
For example, if we require an enhanced concentration to fit data 
for a particular model in OWLS, then haloes measured in this bary- 
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Figure 5. Best fitting halo-model power to the power spectra of the OWLS 
simulations for the DMONLY (black; solid), AGN (purple; long-dashed), 
REF (green; medium-dashed) and DBLIM (red; short-dashed) models up to 
k = 10/zMpc -1 at z = 0.5. These are obtained by fitting both A and rj o 
(equations 14 and 26) to each model at this redshift. In the top panel, we 
show A 2 while in the middle panel we show the ratio of each spectrum to 
the emulator DMONLY case (black crosses in the top panel); one can see 
that the freedom introduced by allowing these parameters to vary is able to 
capture both the up- and down-turn in power that feedback introduces rela¬ 
tive to the dark-matter-only case. Any residual differences for k < 1/zMpc -1 
are due to residual errors in our fitting across a range of cosmologies that 
can be seen in Fig. 2. Our accuracy is best appreciated in the lower panel, in 
which we show the ratio of each halo-model prediction to the corresponding 
simulation. 

onic model should display enhanced concentrations relative to their 
DMONLY counterparts. This approach differs from that presented 
in Semboloni et al. (2013), Fedeli (2014), Fedeli et al. (2014) and 
Mohammed et al. (2014) in that we do not attempt to add accurate 
profiles for the gas and stars into the halo model, but instead look 
for a more empirical modification that is able to match data at the 
level of the power spectrum for k < 10/zMpc -1 . 

Given the above discussion, we might expect that two parame¬ 
ters would suffice: one to capture the increased concentration as gas 
cools in halo cores and one to capture the puffing up of halo profiles 
due to more violent feedback. To fit the baryonic models, we allow 
ourselves to vary the parameter A in the c(M) relation (equation 14) 
and the parameter rjo, where 77 is defined in equation (26) and 170 is 
the first parameter in the full expression for 77 in Table 2, explicitly 

r? = 770 -0.3 ct s (z) • (29) 

All other parameters are fixed to their values in Table 2, and the 
redshift evolution in the second half of the expression for rj is pre¬ 
served from the best fit to the COSMIC EMU nodes. We vary rjo and 
A to best fit the OWLS data from simulations DMONLY, AGN, REF 
and DBLIM. We construct these power spectra by taking ratios of 
the publicly available OWLS baryon models to the DMONLY mod¬ 
els (which produces a smooth curve because the simulations have 



A 

Figure 6. Best matches to the power spectrum from the OWLS simula¬ 
tions found by varying halo structure via A and rjo (equations 14, 26 and 
29) from z = 0 to 1. The contours enclose regions of parameter space that 
match the power spectra with an average error of 2.5 per cent (inner) and 
5 per cent (outer) from k = 0.01 to 10/vMpc and the crosses mark the 
best-fitting point. We show contours for the DMONLY (black; solid), AGN 
(purple; long-dashed) REF (green; medium-dashed) and DBLIM (red; short- 
dashed) cases. These ranges can be used to place a prior on the range of 
770 and A to be explored in a cosmological analysis as they encompass the 
range of behaviour expected from plausible feedback models. The dashed 
line (equation 30) shows a relation between rjo and A that could be used 
to provide a single-parameter fit to all models. The grey cross is the best¬ 
fitting value to all the COSMIC EMU simulations, whereas the black cross is 
the best match to the specific cosmology used in the DMONLY model. 


matched initial conditions) and then multiplying this ratio by the 
COSMIC EMU prediction for the baseline WMAP 3 cosmology used 
for the OWLS simulations. We do this because the OWLS simu¬ 
lations are small in volume and the power spectrum would be too 
noisy to use in its raw form. The best fits to OWLS power spectra 
are shown in Fig. 5 where it can be observed that the freedom per¬ 
mitted by fitting A and rjo allows the power spectrum of HMCODE 
to trace the residual displayed by the OWLS simulations accurately 
over the range of scales shown. Particularly, note that the variation 
is able to reproduce both the up-turn due to gas cooling, enhancing 
clustering around k = lO/iMpc -1 , and the down-turn due to mass 
being expelled from the halo, which can impact the relatively large 
scale of k = 0.3/tMpc . 

In Fig. 6 we show how the goodness-of-fit varies as parame¬ 
ters A and rj are varied for the various feedback recipes. The con¬ 
tours enclose regions of parameter space in which the average error 
is 2.5 per cent (inner) and 5 per cent (outer), where the average 
is taken over all scales between k = 0.01 and lO/iMpc -1 , binned 
logarithmically. One can see that these parameters distinguish well 
between the simulated AGN model and the other models, DBLIM is 
marginally distinguished, but parameters that fitted DMONLY and 
REF best are nearly identical. The distinguishability is directly re¬ 
lated to the magnitude of the effect that each model has on the 
power spectrum (for k < lOfiMpc -1 ), which can be seen in the 
middle panel of Fig. 5. Our best-fitting parameters for each model 
are given in Table 4. The AGN model clearly favours less concen¬ 
trated haloes, which is expected given that AGN blow gas out of the 
central portions of haloes. The range of acceptable parameter com¬ 
binations of A and rjo that are able to fit the OWLS data, shown in 
Fig. 6 could be used to form a prior in future weak-lensing analyses 
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Table 4. Parameter combinations of rjo and A that best fit OWLS data from 
Z = 0 to 1 via the halo-model approach described in the text. These param¬ 
eters are those at the centres of the ellipses in Fig. 6. The OWLS simula¬ 
tions can be matched at the 5 per cent level over the redshift range. That 
the values of t]o and A differ in the case of ‘all COSMIC EMU simulations’ 
compared to DMONLY is because a slightly improved fit is possible in the 
case of dealing with a specific cosmology, which in the case of OWLS 
is the slightly outdated WMAP 3 (£2 m = 0.238, fib = 0.0418, Cg = 0.74, 
n s = 0.951, h = 0.73). 


Model 

no 

A 

All COSMIC EMU simulations 

0.60 

3.13 

DMONLY (WMAP 3 from OWLS) 

0.64 

3.43 

AGN 

0.76 

2.32 

REF 

0.68 

3.91 

DBLIM 

0.70 

3.01 


that aim to constrain or marginalize over baryonic feedback. This 
is acceptable because it is not clear which, if any, of the OWLS 
models is the correct one. The most conservative assumption is that 
the recipes give a range of plausible feedback effects and this is 
what the range we suggest in Fig. 6 encapsulates. Additionally, the 
dashed line shows a relation between rjy and A that could be used 
to provide a single parameter fit to all models: 

7]o = 1.03-0.11A. (30) 

This relation exists because there is some degeneracy between the 
effects of varying rj and those of varying A. Applying this relation 
would mean that the REF model could not be distinguished from 
DMONLY, but all other models could be distinguished. A further 
advantage of the halo-model approach is that it should also cap¬ 
ture any coupling between cosmological parameter variation and 
feedback processes. This effect is ignored in any polynomial- or 
template-based approach to modelling feedback (e.g., Eifleretal. 
2014). 


6 WEAK GRAVITATIONAL LEASING 

The origin of the late-time accelerated expansion of the cosmos is 
uncertain. To test different models it is necessary to measure the 
expansion rate together with the growth of perturbations around 
the present day. The weak gravitational lensing of galaxies has 
emerged as one of the premier tools to probe perturbations, which 
can discriminate between models with similar background expan¬ 
sion rates. 

However, the matter power spectrum is not a directly mea¬ 
surable quantity. With the notable exception of Brown et al. (2003) 
and the 3D lensing of Kitching et al. (2014), the majority of weak- 
lensing analyses have measured real space angular correlations of 
the shears of galaxy images. Shear is typically expressed in terms 
of the spin-2 quantity y, which can be split into tangential (y t ) and 
cross (y x ) components with respect to a pair separation vector. Cor¬ 
relation functions t;± are defined as combinations of two-point cor¬ 
relations of tangential- and cross-shear as a function of angular sep¬ 
aration: 

€± = (rtrt)±(rx7x> • (3i) 

These can be related to the harmonic coefficients, C(£), of the 
spherical Fourier transform of a weighted, projected matter field 


0 [arcminutes] 
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Figure 7. Lensing correlation functions predicted by various methods are 
shown divided by the correlation functions predicted using HMCODE with 
no k-cut imposed. §+ ( red) is shown in the upper panel and (or¬ 
ange) in the lower panel. Source galaxies are all taken to be fixed at 
z s = 0.7, the effective median lensing redshift of CFHTLenS. We show the 
cases of taking the HMCODE matter power spectrum to be sharply cut at 
k— 10h (short-dashed) and 50/zMpc 1 (long-dashed) and correlation func¬ 
tions from HALOFIT (dot-dashed) and COSMIC EMU (solid; where A 2 = 0 
for k > 10/iMpc -1 ). For the range of angular scales shown, all models agree 
for §+ at the ~ 3 per cent level, a k-cut at k — 10/iMpc -1 only starts to im¬ 
pact upon £+ at the per cent level for 6 < 0.03° and a cut at 50/zMpc -1 has 
almost no impact for angles shown; the agreement between HMCODE and 
COSMIC EMU is ~ 2 per cent across the range of angles. For HMCODE 
and COSMIC EMU agree to ~ 2 per cent until the k- cut becomes important, 
around 6 = 0.2°. HALOFIT is discrepant at the 4 per cent level for 6 < 0.5°. 
For , the impact of k-cuts is felt at larger angular scales than for - a 
cut at 50/zMpc -1 makes its impact at the per cent level around 6 = 0.02°. 


(Kaiser 1992), 

£t(0) = T rc(l)J±mt df , (32) 

2n J o 

with harmonic wavenumber l. Here 7± is the zeroth (£+) and fourth 
(£_) Bessel functions respectively. In the Limber (1954) approxi¬ 
mation the C(£) are then related to integrals over all scales of the 
matter power spectrum (Kaiser 1992): 

C{£) = r = *//*(®).z(®)) d “ . (33) 

Jo a z (a>) 


where ffl(z) is the comoving distance to redshift z, defined via the 
metric convention of Bartelmann & Schneider (2001) and g(co) is 
the lensing efficiency function: 


S(®) 


3H° 5 Q m /•<** frJ J K (co'-(o) 

-&-L 


dco 1 . 


(34) 


Here /jf (tn) is the comoving angular-diameter distance, 0 )h is the 
comoving distance to the horizon and ra(co) is the normalized dis¬ 
tribution of source galaxies. Finally, P(k) is the power spectrum of 
matter fluctuations, defined in equation (4), and is exactly what HM¬ 
CODE provides. Thus, to make contact with lensing observables, it 
is necessary to investigate how the accuracy of the predictions of 
HMCODE at the level of the matter power spectrum translates into 
accuracy for 

Several methods are commonly used to deal with the fact that 
2D lensing correlation functions depend on integrals over all k of 
the matter spectrum. Either fitting formulae are extrapolated be¬ 
yond the regime to which they were fitted, in the hope that they 
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Figure 8. As Fig. 7 but for the C(l ) coefficients. Each model shown is 
divided by the HMCODE prediction with no k-cut imposed. We see that a cut 
in A 2 at k = 10/? Mpc 1 induces percent level deviations around £ = 10 3 and 
a cut at 50//Mpc 1 increases this to £ — 10 4 . HMCODE and COSMIC EMU 
agree at the per cent level until the finite /.'-range of COSMIC EMU becomes 
important. HALOFIT disagrees by as much as 5 per cent around i = 10 3 and 
this disagreement increases to 7 per cent for higher harmonics. 

still provide the required accuracy, or P(k ) is set to zero for re¬ 
gions where the power is unknown, a so-called k- cut. As previ¬ 
ously stated, one of the benefits of our halo-model approach is that 
we have more reason to trust the predictions of the halo model 
at small scales, due to the large amount of theoretical input that 
goes into the model. In Fig. 7 we show the lensing t,± correla¬ 
tion functions as predicted by integrating over P(k) from HMCODE, 
HALOFIT and from COSMIC EMU (where we set A 2 (k) = 0 for 
k > lOfiMpc -1 ). These are computed by taking the source galaxies 
to all be fixed at z s = 0.7, the effective median redshift for lens¬ 
ing of CFHTLenS (Fleymans et al. 2012), and using the best-fitting 
cosmology from the Planck Collaboration XIII (2015). The three 
models for agree at the ~ 3 per cent level for at all angular 
scales shown with predictions from HMCODE and COSMIC EMU 
agreeing at the per cent level for all angles shown. We also show 
HMCODE <!;± predictions for k-cuts at 10 and 50/7 Mpc -1 . If the the¬ 
ory were perfect, a cut-off in A 2 at k = 1 0/7 Mpc -1 provides per cent 
level accuracy only for 6 > 0.03°; if the power spectrum is cut at 
50/7 Mpc -1 , then the predicted t,+ does not deviate from per cent 
level agreement for all angles shown. For if no k-cut is imposed, 
the three models agree at the per cent level only for 8 > 1°, after 
which predictions from Takahashi et al. (2012) deviate by as much 
as 6 per cent at 8 ~ 0.01°. The effect of the finite resolution of COS¬ 
MIC EMU becomes important at the per cent level at 8 = 0.2° and 
the same deviation can be seen in the HMCODE prediction when it 
is cut at k = 1 Ofi Mpc -1 . Assuming perfect knowledge of the theory 
this accuracy extends to 8 > 0.02° if the cut is taken at 50/7 Mpc -1 . 
We note that if we extend the matter spectrum from COSMIC EMU 
using a power law for k > lOftMpc -1 the predictions agree with 
those of HMCODE to 2 per cent for both %±; this is probably due to 
the fact that the halo-model prediction at k > 10/? Mpc -1 involves 
an integral over many quantities that are accurately power laws, 
resulting in a close to power-law power spectrum. 

In Fig. 8 we compare models at the level of their C(£) pre¬ 
dictions. When the finite k-range of COSMIC EMU is unimportant, 
it agrees with the HMCODE prediction to one per cent. Discrepan¬ 
cies with HALOFIT at the 5 per cent level arise from t = 500 with 
maximum deviations of 7 per cent for l > 10 4 . A cut in power at 
k = 10/;Mpc -1 impacts upon the C{tj at the per cent level around 
l = 10 3 and a cut at 50/iMpc -1 at t = 10 4 . 

Many authors have investigated how baryonic processes af¬ 
fect weak-lensing observables. Early work (White & Vale 2004; 
Zhan & Knox 2004) used the halo model to estimate how the mat¬ 
ter power would be altered by including gas cooling in haloes 
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Figure 9. (upper panel) and c (lower panel) correlation functions pre¬ 
dicted using HMCODE, with a width showing the spread that is obtained 
from different feedback models. In each case, source galaxies are taken to 
be fixed at z s — 0.7, approximately the effective median redshift for lensing 
for CFHTLenS, and we show the correlation functions predicted using the 
best fitting Planck cosmology (upper curve) and the best fitting CFHTLenS 
(Heymans et al. 2013) cosmology (lower curve). For comparison, we also 
show the measured t;± from the CFHTLenS survey; it can be seen that feed¬ 
back fails to alleviate the tension between CFHTLenS and Planck data. One 
can see that an ignorance of the details of feedback affects g to much larger 
angular scales than g t. a consequence of it probing more non-linear regions 
of the matter distribution. Baryonic feedback has an impact at the greater 
than per cent level for 0 <0.1° for (L and 9 <2° for c, . In all cases, the 
effects of baryonic feedback are small relative to the errors in current data. 

and hot, diffuse intra-cluster gas, respectively. More recent work 
has used hydrodynamic simulations with feedback recipes (e.g., 
Jing et al. 2006; Semboloni et al. 2011; Casarini et al. 2012) to 
compare weak-lensing observables to the case when no baryonic 
feedback is included. The results are that for £ > 1000, the C(£) are 
altered at the per cent level with the alterations increasing with (, 
but the details are strongly dependent on the feedback implementa¬ 
tion. 

In Fig. 9 we show the range of possible correlation function 
predictions given by our power spectrum fits to the baryonic feed¬ 
back models, where the region enclosed by the curves is the region 
that our fits to the OWLS feedback models occupy (the centres of 
the ellipses in Fig. 6). We generate these using HMCODE predic¬ 
tions with parameters A and rjo taken from the centres of the el¬ 
lipses in Fig. 6. We also show data from the CFHTLenS analysis of 
Kilbinger et al. (2013) so that ignorance of the details of feedback 
can be compared to the current errors in data. For current data we 
see that the effect of feedback is small compared to the errors, but 
data that will be available in the near future will increase in accu¬ 
racy and feedback processes will have to be accounted for. Fig. 9 
also shows that baryonic feedback does little to alleviate the tension 
between the best-fitting CFHTLenS cosmology and that of Planck. 
For our case of sources fixed at z s = 0.7 baryonic feedback only has 
an effect at the greater than per cent level for 8 < 0.1° for and 
8 < 1° for In a forthcoming paper (Joudaki et al. in prep) con- 
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straints on A and 7 ]q, together with cosmological constraints when 
these parameters are marginalized over, will be presented using the 
CFHTLenS together with that from RCSLenS. 

Alternative approaches have been investigated to model 
the impact of feedback on weak-lensing observables, all of 
which use data from the OWLS hydrodynamic simulations: 
Mohammed & Seljak (2014) model the OWLS data by refitting 
the coefficients from their power-series expansion of the one-halo 
term and advocate marginalizing over these coefficients to immu¬ 
nize against biases due to feedback. Harnois-Deraps et al. (2015) 
construct polynomial fits to the ratio of power spectra from feed¬ 
back models to the DMONLY model; again the coefficients of these 
polynomials could feasibly be constrained by data. However, to fit 
each model over the scale redshift range required 15 coefficients, 
compared to only 2 in our approach. With a fixed (WMAP 9) cos¬ 
mology, Harnois-Deraps et al. (2015) find a preference for feed¬ 
back in the CFHTLenS data. MacCrann et al. (2015) reanalysed the 
CFHTLenS survey but adding a single parameter that governs the 
amplitude of AGN feedback, which was taken to be given by the 
ratio of power from the AGN to DMONLY simulations. They find 
only a weak preference for feedback, but find that AGN feedback 
is insufficient to resolve tension between the CFHTLenS data and 
that of the Planck satellite. An identical conclusion is reached by 
Kitching et al. (2014) using a similar method. Eifleretal. (2014) 
propose using a principal component analysis method, by which 
components of the power spectrum that are most affected by feed¬ 
back are removed. Assuming that feedback is independent of cos¬ 
mology, they show that the data could be fitted with as few as 4 
components removed. Our method may be preferable to this as it 
potentially allows one to capture the coupling between feedback 
and cosmological parameters, via their effects on the halo profiles. 


7 SUMMARY AND DISCUSSION 

We have shown that the halo model can be optimized so that it ac¬ 
curately reproduces power spectra measured from dark-matter N- 
body simulations across a range of cosmological models for k < 
10/7Mpc~’ and z < 2, provided we are willing to introduce a num¬ 
ber of empirical modifications to its ingredients. We achieved this 
by calibrating our model to the ‘node’ simulations of COSMIC EMU 
of Heitmann et al. (2014). Our success reflects the fact that the halo 
model is built on well-posed theoretical ingredients, which natu¬ 
rally adapt to changes in cosmology in a sensible fashion. Our fits 
are accurate at the 5 per cent level, which represents an improve¬ 
ment over the currently used HALOFIT model of Takahashi et al. 
(2012). COSMIC EMU itself is quoted to be accurate to 5 per cent 
at k = lO/iMpc -1 and so our fit is as good as possible at the cur¬ 
rent level of ignorance at this scale. Even in the dark-matter-only 
case, our accuracy statement comes with the caveat that we have 
only tested a limited range of plausibly interesting cosmologies. In 
particular, we have concentrated on the parameter cube of the COS¬ 
MIC EMU of Heitmann et al. (2014), which contains cosmological 
parameters within the 5 a region of the Planck Collaboration XIII 
(2015) results for the standard cosmological model. 

The advantage of the halo-model approach is that it can be 
simply expanded beyond the parameter cube of COSMIC EMU. We 
demonstrated that our model is able to produce reasonable power 
for k > lOfi Mpc~* and it can also produce spectra for z > 4, 
greater than allowed by COSMIC EMU. Given the large amount 
of tested theoretical input that goes into the halo-model calcula¬ 
tion, we expect that our model should produce sensible spectra for 


higher redshifts and smaller scales. For small deviations from the 
standard cosmological paradigm, such as dark energy with time- 
varying equation of state, or for models with small amounts of cur¬ 
vature, we also expect our answers to be accurate. If one were in¬ 
terested in the power spectrum of radically different models, such 
as very curved models or those qualitatively different linear power- 
spectrum shapes, we would advise caution, as our model has not 
been tested to these extremes. 

Our approach differs from that of some authors (e.g., 
Cooray & Hu 2001; Giocoli et al. 2010) who attempt to improve 
the basic halo model by adding effects such as halo substructure 
and a scatter in halo properties at a given mass. It may be that 
adding in these ingredients would have reduced the number or 
magnitude of our fitted parameters, but we make no attempt to 
quantify this. Our approach also contrasts with that employed by 
Mohammed & Seljak (2014) and Seljak & Vlah (2015) who ad¬ 
vocate replacing a theoretically motivated one-halo term with a 
power-series expansion and fitting this expansion term-by-term. 
Combining this approach with perturbation theory produced excel¬ 
lent results in the quasi-linear ( k ~ 0.3/iMpc -1 ) regime. However, 
it is not obvious how to extend their predictions to smaller scales 
than the smallest constrained by these authors (k ~ lfiMpc - ), as 
their empirical power series has no physical requirement for sensi¬ 
ble behaviour at smaller scales. Given the necessity of these smaller 
scales in producing usable weak-lensing predictions, we therefore 
prefer our approach. A future avenue of fruitful research may be to 
combine the power of our approach of a physically motivated fit to 
the one-halo term with the perturbation-theory-inspired two-halo 
term advocated in these works. Particularly this may improve ac¬ 
curacy in the quasi-linear regime and would certainly help in mod¬ 
elling around the BAO scale, where our 5 per cent accuracy deviates 
from the per cent level accuracy of COSMIC EMU. 

Baryonic feedback has a large impact on the power at small 
scales. We demonstrated that the halo model is able to capture the 
influence of baryonic physics via only two redshift-independent pa¬ 
rameters that govern halo internal structure. We also showed that 
this can be reduced to one parameter, at the loss of a small amount 
of discriminating power. Using these two parameters, we were able 
to fit feedback recipes considered in the OWLS simulations at the 5 
per cent level for j < 1 and k < 10/jMpc~’. Because these param¬ 
eters are firmly rooted in the halo-model apparatus, their effects 
are restricted to scales at which haloes affect the power spectrum. 
This is not guaranteed by other approaches, such as polynomial 
fits, which may produce unphysical effects. We also suggested that 
our approach is more likely to reproduce the correct coupling be¬ 
tween baryonic feedback and cosmology, because it is rooted in 
halo properties. It is not obvious how to account for the cosmology 
dependence of baryonic feedback using existing fitting formulae or 
COSMIC EMU. 

Finally, we showed how the power-spectrum predictions of 
HMCODE translate into the lensing C(£) and correlation func¬ 
tions that are measured in a standard lensing analysis. For the dark- 
matter-only case, we showed that HMCODE agrees with COSMIC 
EMU at the per cent level for E, + and C(£) and 2 per cent for for 
all scales where the finite Urange of COSMIC EMU is unimportant. 
We suggest that HMCODE provides a reasonable way of extrapolat¬ 
ing COSMIC EMU to the smaller scales that are necessary to pro¬ 
duce weak-lensing predictions. Considering baryonic feedback, we 
showed how our matter power spectra translate into lensing pre¬ 
dictions and showed that ignorance of the details of feedback is 
smaller than uncertainty on current data. For a survey similar to 
CFHTLenS, baryonic feedback impacts on at the per cent level 
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for 6 < 0.1° and on for 6 < 2°. In future tensing analyses (e.g., 
Joudaki et al., in preparation) we advocate marginalising over the 
halo parameters that we used to fit to the OWLS feedback models 
in order to produce unbiased cosmological constraints; our range 
of fitted values for these parameters may be used as a prior for this 
purpose. Alternatively, one might accept a given cosmology and 
then use the best-fitting baryon parameters as a means of learning 
about baryonic feedback. 

In summary, and given our accuracy, we suggest that COS¬ 
MIC EMU be used if one is interested in the non-linear, gravity- 
only induced power spectrum for k < lO/iMpc -1 , and the required 
model is within the COSMIC EMU parameter cube. However, if one 
is interested in departures from the COSMIC EMU parameter space, 
accounting for the effect of baryonic feedback physics, or produc¬ 
ing accurate lensing observables via a reasonable extrapolation, we 
then advocate HMCODE. Although we focused on weak lensing in 
this paper, we stress that HMCODE is useful for any application that 
currently uses HALOFIT. 

This last point emphasizes the potential of the approach de¬ 
scribed in this paper. The halo model can readily be extended to 
take account of new physical processes and changes in the cos¬ 
mological paradigm. One example for further work would be an 
application of our method to cover modified gravity models (e.g., 
Schmidt etal. 2010; Lombriser et al. 2014; Barreira et al. 2014) 
where revised growth rates, collapse thresholds and internal halo 
structures can be predicted in part on analytic grounds, and where 
there is a growing effort on detailed simulations. Another such ex¬ 
ample would be to look at the impact of massive neutrinos (e.g., 
Massara et al. 2014) on the matter spectrum. In such cases, being 
able to produce accurate power spectra will be important in or¬ 
der to distinguish standard and non-standard cosmological models. 
Moreover, exploration of a large parameter space of models will in¬ 
evitably be necessary, and there will therefore be a strong motiva¬ 
tion to explore rapid means of generating non-linear power spectra. 

The code developed as part of this paper is available at 
https://github.com/alexander-mead/hmcode or at request 
from the author. It is able to produce the matter power spectrum at 
200 k values for 16 different z values in ~ 0.5s on a single core of 
a standard desktop computer. 
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APPENDIX A: POWER SPECTRUM RESPONSE 

In Fig. Al we show the response of the matter power spectrum to changes in 
cosmological parameters for the linear spectrum, and the non-linear spec¬ 
trum predicted by both HMCODE and COSMIC EMU. This is shown for the 
cosmological parameters that COSMIC EMU allows one to vary (co m , (0 b, 
n s , w, as h ) over the range that COSMIC EMU allows (see Table 1) via 
the ratio of the power as the cosmological parameter is varied compared 
to the power for the ‘fiducial’ cosmology, given in Table 1. This quantity 
says nothing about the absolute accuracy of the HMCODE predictions, but 
allows us to assess if the response of the power induced by changes in cos¬ 
mological parameters is accurate. Therefore, we also compute the ratio of 
the response of the HMCODE and COSMIC EMU predictions (right-hand col¬ 
umn). We see that, in general, the matter power spectrum response pre¬ 
dicted by HMCODE is in excellent agreement with the simulations of COS¬ 
MIC EMU, which is due to the large amount of well-tested theory in the 
halo model. The most obvious small discrepancies arise around the ‘bump’ 
at k ~ I /2 Mpc 1 produced by variations in as and the degree of enhanced 
power at k > l/?Mpc _l when w is increased. In each case, the general trend 
is reproduced well by the halo model but the exact magnitude of the re¬ 
sponse is not predicted quite correctly. 

The accuracy of the predicted halo-model power spectrum depends 
on the ingredients used. We used Fig. Al to inform the ingredients used for 
HMCODE in order to produce the correct power-spectrum response, before 
we fitted halo parameters to make accurate A 2 (k) predictions. This was par¬ 
ticularly important in the case of the response to variations in w, which enter 
our model partly through the concentration-mass relation of Bullock et al. 
(2001) and partly via the Dolag et al. (2004) correction to a c(M) relation. 
The vast majority of concentration-mass relations available in the literature 
produce either no, or too little, response in the matter spectrum to changes 
in w when the linear theory A 2 (£) is held constant. 


APPENDIX B: THE HALO MODEL CALCULATION 


In this appendix we discuss some of the practicalities of our approach of cal¬ 
culating the halo-model power spectrum. This appendix serves as a manual 
to the specific form of the halo-model calculation that our public FORTRAN 
90 code performs. The halo-model code used in this work is available at 
https://github.com/alexander-mead/hmcode. 

HMCODE is written such that by default it uses the Eisenstein & Hu 
(1998) approximation to the linear spectrum, which can be used if accuracy 
at linear scales is not demanded, and then converts this to a non-linear spec¬ 
trum. If necessary there is the option to read in a tabulated linear spectrum, 
k versus P{k), as input (e.g., from CAMB). Additionally, the cosmological 
parameters need to be specified because they are used in the fitting func¬ 
tions and in calculations of the growth function. In all cases the input linear 
spectra are renormalized to give the desired as. 

For our calculation of the growth function, we explicitly integrate the 
approximate expression for the logarithmic growth rate given in Linder 
(2005) and derived in Linder & Cahn (2007), 


ding _ y , , 

11 — ^m(z) • 

alna 


(Bl) 


where g is the growth factor normalized to be 1 today, y — 0.55 if w — — 1 
and y= 0.55 + 0.02(l+w) if w < — 1 and y= 0.55 + 0.05(1 +w) if w > —1. 
This fitting formula and subsequent integration to find the growth factor is 
valid at the sub-per cent level even if w deviates significantly from — 1. 
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Figure Al. The ratio of matter power at z = 0 as cosmological parameters are varied when compared to a ‘fiducial’ cosmological model of the COSMIC EMU 
parameter space. We show comparisons of linear theory (left column), HMCODE (second column) and the COSMIC EMU prediction (third column). Comparing 
the two central columns, one can see that in general HMCODE is able to reproduce the trends in power-spectrum response to cosmological parameter variation 
accurately. The right-hand column shows a residual of the HMCODE power divided by that of the emulator (i.e. each curve in the second column divided by the 
corresponding curve in the third column), which draws attention to the response of the power spectrum that is replicated best and least well by HMCODE. In 
each case, the cosmological parameter in question is varied over the range that COSMIC EMU covers, while all other parameters are kept fixed, with the highest 
value of the cosmological parameter shown in pink and the lowest in blue; the range of each parameter is given in Table 1. This plot is slightly different from 
similar plots in Heitmann et al. (2014) because we use k/h rather than k on the x-axis. 
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B1 Two-halo term 

The two halo term we use is: 

A?h« = [l ~/tanh 2 (tov/V?)] A 2 in (k) , 

with / = 0.188 x <jg' 29 (z). The parameter <7 V is calculated via 


u 


■dk, 


(B2) 


(B3) 


3 Jo k 3 

and we do this by transforming the [0,«>] interval in k to t G [0,1] using the 
transformation 1 + k = 1 ft. 

The Ov integral (and that of <7 (R) in equation B5) requires knowledge 
of A 2 at very small or very large values of k. It is impractical to provide a 
tabulated linear spectrum over such a large range of scales and so we inter¬ 
polate beyond the boundaries of the input linear spectrum using the scaling 
A 2 (£) oc £ 3 +”s a t small k, and the approximate scaling A 2 (k) °c £" s_1 \n 2 (k) 
at high k. Extremely accurate values of the power at each of these asymp¬ 
totes are not necessary, because they contribute quite negligibly to the inte¬ 
gral, but it is necessary for the linear power not to be badly wrong, or set to 
zero unphysically. 


virial radius, related to the halo mass via M — 4^r 2 A v p/3 and A v = 
418 x a-° 352 (z). The halo concentration is calculated using the prescrip¬ 
tion of Bullock et al. (2001), augmented by that of Dolag et al. (2004): 

C (M,;)=A 1 i + ZfgDE / Z ^°° ) , (B8) 

1+2 8a(z ^ 00 ) 

where A — 3.13. The growth factor correction only applies if w ^ — 1 and 
is the asymptotic ratio of growth factors. The formation redshift is then 
calculated via 

s(Zf). 


s(z) 


o{fM,z) = 8 C 


(B9) 


where / = 0.01. We numerically invert equation (B9) to find if for a fixed 
M. If Zf < z, then we set c = A. The mass function we use in equation (B4) 
is that of Sheth & Tormen (1999): 

1 (BIO) 


/(v)=A 


1 - 


-av 2 /2 


(av 2 )''J 

where a = 0.707, p = 0.3, A = 0.2162 and V = 8 C /<J(M) with 8 C = 1.59- 
0.0314 In og(z). 


B2 One-halo term 

The aim is to numerically evaluate the integral in equation (8). To do this 
we find it convenient to convert from an integral over M to one over V = 
8 c /a(M): 

, 3 


Ahm =[i(4) i 

x [ M(v)W 2 (v^k,M)f(v) dv . 
Jo 


(B4) 


Here k * = 0.584 a/ 1 (z) and damps the one-halo term to prevent one-halo 
power from rising above linear at the largest scales, r] — 0.603 — 0.3 <78 ( 2 ) 
and bloats haloes while they maintain a constant virial radius. We integrate 
over a finite range in v that captures all haloes necessary to produce a con¬ 
vergent power spectrum at the scales we investigate. In testing we found that 
v G [0.3,5] was sufficient to produce convergent results to k = 100/iMpc -1 
and v G [0.1,5] if one required power out to k = 10 4 /iMpc _1 . Convergence 
at higher wave numbers requires the minimum value of v to be reduced, 
because low-mass haloes contribute to the power only at small scales. 

The halo-model integral in equation (B4) requires knowledge of r v , 
M, c and a, all as a function of halo mass. In practice, we tabulate these 
over the finite range in v as an ‘initialization’ step to the calculation. We 
optimized the number of points in v in order to produce convergent results 
up to k — 100/zMpc -1 . o(R) is computed via 


e 2 (R)= f A 2 (k)T 2 (kR)d\nk . 
Jo 

where 

3 

T(x) — -^-(sinx — xcosx) . 


(B5) 


(B6) 


We convert the [0,°o] k range of the integral in equation (B5) to a finite 
range in t using 1 -\-k — \/t with t G [0,1]. Since this integral is relatively 
expensive to compute, we generate a look-up table for <7 (R) with R values 
logarithmically spaced between 10 -4 and 10 3 /z -1 Mpc. These radii corre¬ 
spond to all haloes of practical interest at low z • If values of <7 (R) are re¬ 
quired outside this range, then we interpolate beyond the table boundaries. 
This works extremely well because <7 (R) is a very smooth function in log 
space with the asymptotes approximating power laws to a very high degree 
of accuracy for standard cosmological spectra. The radial scale R is related 
to the mass scale via M = 4nR 2 p/3 . 

The halo window function in equation (B4), W(k,M), has an analyti¬ 
cal form for an NFW profile (Cooray & Sheth 2002): 

W(k,M)f(c) = [Ci(£ s (l +c)) — Ci(fc s )]cos(fc s ) 

sin(cL) (B7) 

+ [Si(*.(l +c)) - Sift)] Sinft) - , 

where /(c) = ln(l +c) — c/(l +c), Si(x) and Ci(x) are the sine and co¬ 
sine integrals, k s = kr v /c, c is the halo concentration and r v is the halo 


B3 Full power 

The final expression for the halo-model power spectrum is then 
A 2 (i) = [(A' 2 2 1 )“ + (A' 2 H )“] 1 /“, 
where a = 2.93 x 1.77" cff and 
din O' 2 (R) 

( 7=1 


(Bll) 


3+«eff= -- 


dlnR 


(B12) 


We calculate the derivative numerically from our (7 (R) look-up table. 
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